Tutorial Part VI — Population analysis in R

Introduction

This is a brief tutorial about constructing and analyzing relatively simple Cormack-Jolly-Seber capture-mark-recapture (CMR) models in R. For this we use the R package “marked” (Laake et al. 2013).

First we get an overview of the data, then we develop different models and finally we analyse and visualize the results.

Setup

Load (and install) the following packages:

package.list=c("here",
               "marked",
               "skimr",
               "tidyverse",
               "sf",
               "tmap"
               )

for (package in package.list) {
  if (!require(package, character.only=T, quietly=T)) {
    install.packages(package)
    library(package, character.only=T)
  }
}

set the theme for some ggplots:

theme_set(
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    panel.border = element_rect(
      colour = "black",
      fill = NA,
      linewidth = 1
    ),
    axis.ticks = element_line(colour = "black", linewidth = 1),
    axis.ticks.length = unit(0.2, "cm")
  )
)

Load and manipulate data

Here, we want to directly set some columns as factors for all subsequent analyses:

starlings <-
  read_csv(here("data", "data_berlin", "animal_data", "starlings.csv")) %>%
  mutate(
    sex = as.factor(sex),
    habitat = as.factor(habitat),
    infected = as.factor(infected)
  ) %>%
  rename(ID = "...1")

# Short overview
skim(starlings)
Data summary
Name starlings
Number of rows 605
Number of columns 5
_______________________
Column type frequency:
character 1
factor 3
numeric 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
ch 0 1 7 7 0 32 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
sex 0 1.00 FALSE 2 Fem: 323, Mal: 282
habitat 0 1.00 FALSE 2 urb: 311, rur: 294
infected 202 0.67 FALSE 2 Yes: 214, No: 189

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
ID 0 1 303 174.79 1 152 303 454 605 ▇▇▇▇▇ 
set.seed(123)

Each row in this table represents one individual. We have 5 columns in this dataset: ID, the individual identifier, ch, the capture history (0 = undetected, 1 = detected),the sex of the captured individuals, the habitat type the birds were captured, and health status, i.e. if birds were infected

Explore capture locations

#load environmental data on capturing
locations <-
  read_csv(here(
    "data",
    "data_berlin",
    "animal_data",
    "starlings_capture_location_data.csv"
  ))

# Short overview
skim(locations)
Data summary
Name locations
Number of rows 2
Number of columns 10
_______________________
Column type frequency:
character 1
numeric 9
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
Habitat 0 1 5 5 0 2 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Lon 0 1 13.57 0.19 13.44 13.51 13.57 13.64 13.71 ▇▁▁▁▇
Lat 0 1 52.44 0.05 52.40 52.42 52.44 52.45 52.47 ▇▁▁▁▇
temperature_2016 0 1 30.25 1.06 29.50 29.88 30.25 30.62 31.00 ▇▁▁▁▇
temperature_2017 0 1 39.00 5.66 35.00 37.00 39.00 41.00 43.00 ▇▁▁▁▇
temperature_2018 0 1 33.50 3.54 31.00 32.25 33.50 34.75 36.00 ▇▁▁▁▇
temperature_2019 0 1 28.50 2.12 27.00 27.75 28.50 29.25 30.00 ▇▁▁▁▇
temperature_2020 0 1 28.50 0.71 28.00 28.25 28.50 28.75 29.00 ▇▁▁▁▇
temperature_2021 0 1 30.00 1.41 29.00 29.50 30.00 30.50 31.00 ▇▁▁▁▇
temperature_2022 0 1 29.50 1.41 28.50 29.00 29.50 30.00 30.50 ▇▁▁▁▇ 
loc_sf <- st_as_sf(locations, coords = c("Lon", "Lat"), crs = 4326)


tmap_mode("view")
tm_shape(loc_sf) + tm_dots(col = "Habitat", size = 2)

Plot the capturing data

# Capture histories
starlings %>%
  ggplot(mapping = aes(x = ch)) +
  geom_bar(colour = "#262674", fill = "#262674") +
  scale_y_continuous(
    limits = c(0, 40),
    breaks = seq(0, 40, by = 10),
    expand = c(.001, .001)
  ) +
  labs(x = "capture histories") +
  theme(axis.text.x = element_text(angle = 45))

# Sex
starlings %>%
  ggplot(mapping = aes(x = sex, fill = sex)) +
  geom_bar() +
  scale_fill_manual(values = c(Female = "salmon", Male = "lightblue"))

# Sex and capture histories
starlings %>%
  ggplot(mapping = aes(x = ch, fill = sex)) +
  geom_bar(position = "stack") +
  scale_y_continuous(
    limits = c(0, 40),
    breaks = seq(0, 40, by = 10),
    expand = c(.001, .001)
  ) +
  labs(x = "capture histories", fill = "sex") +
  theme(axis.text.x = element_text(angle = 45)) +
  scale_fill_manual(values = c(Female = "salmon", Male = "lightblue"))

Plot environmental data

locations_long <- locations %>%
  pivot_longer(names_to = "year",
               values_to = "temp",
               cols = starts_with("Temp"))

ggplot(locations_long, aes(x = year, y = temp, group = Habitat)) +
  geom_line(aes(col = Habitat), linewidth = 2) +
  scale_colour_manual(values = c(urban = "darkgoldenrod", rural = "darkblue"))

Basic Cormack-Jolly-Seber models

The default model uses constant detection and survival probabilities and constant time intervals between captures. We fit the model and have a look at the results. Please note: The estimates are on a logit scale

#models require dataframes, tibbles produce errors
starlings <- starlings %>%
  as.data.frame()

cjs_m1 <- crm(data = starlings)
cjs_m1
## 
## crm Model Summary
## 
## Npar :  2
## -2lnL:  1291.825
## AIC  :  1295.825
## 
## Beta
##                   Estimate
## Phi.(Intercept) 0.06757844
## p.(Intercept)   2.22899510

Npar = number of parameters -2lnL = AIC = Akaikes Information Criterion Phi = Survival probability between capture events p = Capturing probability per capture event

Estimate parameters on real scale

Write a function to assess results more intuitively on the real scale instead of the logit-scale. We can just paste our results from the previous model output

inverse_logit <- function(x) {
  exp(x) / (1 + exp(x))
  
}

#e.g.
inverse_logit(0.06757844) # this is our survival probability
## [1] 0.5168882
inverse_logit(2.22899510) # this is our survival probability
## [1] 0.9028232
#  alternative way to write this

inverse_logit(cjs_m1$results$beta$Phi)
## (Intercept) 
##   0.5168882
inverse_logit(cjs_m1$results$beta$p)
## (Intercept) 
##   0.9028232

Refit the model with confidence intervals

cjs_m1_2 <- cjs.hessian(cjs_m1)
cjs_m1_2
## 
## crm Model Summary
## 
## Npar :  2
## -2lnL:  1291.825
## AIC  :  1295.825
## 
## Beta
##                   Estimate         se         lcl       ucl
## Phi.(Intercept) 0.06757844 0.07249005 -0.07450207 0.2096589
## p.(Intercept)   2.22899510 0.24951951  1.73993685 2.7180533

Predict

Now we transform the logit estimates to real estimates by using the inverse logit

predict(object = cjs_m1_2,
        newdata = data.frame(sex = c("Female", "Male")),
        se = TRUE)
## $Phi
##   occ  estimate         se       lcl       ucl
## 1   1 0.5168882 0.01810184 0.4813831 0.5522236
## 
## $p
##   occ  estimate         se      lcl       ucl
## 1   2 0.9028232 0.02189121 0.850679 0.9380836

What does that mean? The survival probability between capture events was 0.52. The probability to capture an individual during a capture event was 0.90.

As mentioned before, we assume constant time intervals between capturing events. In reality, this is often not the case, e.g. capturing events have to be postponed due to weather conditions. Now we build a model where we take this into account.

CJS with irregular sampling intervals

cjs_m1_time <- crm(data = starlings,
                   # time.interval vector is the interval between each capture event.
                   # for 7 capture events, there are 6 intervals
                   time.intervals = c(1, 1, 1, 1.6, 2, 3))
predict(object = cjs_m1_time)
## $Phi
##   occ estimate
## 1   6 0.688108
## 
## $p
##   occ  estimate
## 1   7 0.8858232

Now we have a higher survival probability of 0.69 between capturing events and a lower capturing probability per capturing event of 0.89

Fit models for model selection

Normally, there are different hypotheses about the factors that influence the survival probability of individuals. These can be characteristics of the individuals, such as sex, age or weight, or characteristics of the environment of the individuals, such as sealing, temperature or the number of nesting opportunities. For this purpose, different models are built that take the respective factors (explanatory or predictor variables) into account.

For a flexible approach we process the data, create design data to add variables that depend on time and fit the different models.

Group data

We process the data and set the grouping variables

m2_proc <- process.data(data = starlings,groups = c("sex","habitat","infected"))

Create the design data

m2_design <- make.design.data(data = m2_proc)

Model formulas

Build formulas for each parameter

phi_dot <- list(formula = ~ 1)
phi_sex <- list(formula = ~ sex)
p_sex <- list(formula = ~ sex)

Fit a model based on the design data with constant survival and different detection probabilities between sexes.

cjs_m2 <- crm(
  m2_proc,
  m2_design,
  model.parameters = list(Phi = phi_dot,
                          p = p_sex),
  accumulate = FALSE
)
## 
 Number of evaluations:  100  -2lnl: 1289.126943
cjs_m2
## 
## crm Model Summary
## 
## Npar :  3
## -2lnL:  1289.117
## AIC  :  1295.117
## 
## Beta
##                  Estimate
## Phi.(Intercept) 0.0793894
## p.(Intercept)   1.7877066
## p.sexMale       0.7906268
predict(object = cjs_m2)
## $Phi
##   occ  estimate
## 1   6 0.5198369
## 
## $p
##      sex occ  estimate
## 1 Female   7 0.8566459
## 2   Male   7 0.9294541

Fit another model based on the design data with different survival and detection probabilities between sexes

cjs_m3 <- crm(
  m2_proc,
  m2_design,
  model.parameters = list(Phi = phi_sex,
                          p = p_sex),
  accumulate = FALSE
)
## 
 Number of evaluations:  100  -2lnl: 1283.680512
 Number of evaluations:  200  -2lnl: 1283.680512
cjs_m3
## 
## crm Model Summary
## 
## Npar :  4
## -2lnL:  1283.681
## AIC  :  1291.681
## 
## Beta
##                    Estimate
## Phi.(Intercept) -0.09718598
## Phi.sexMale      0.34094476
## p.(Intercept)    1.98842947
## p.sexMale        0.46984162
predict(object = cjs_m3)
## $Phi
##      sex occ  estimate
## 1 Female   6 0.4757226
## 2   Male   6 0.5606397
## 
## $p
##      sex occ  estimate
## 1 Female   7 0.8795769
## 2   Male   7 0.9211642

Fit another model based on the design data with time-dependent survival and detection probability

cjs_m4 <- crm(
  m2_proc,
  m2_design,
  model.parameters = list(Phi = list(formula =  ~ time),
                          p = list(formula =  ~ time)),
  accumulate = FALSE
)
## 
 Number of evaluations:  100  -2lnl: 1244.670707
 Number of evaluations:  200  -2lnl: 1239.725494
 Number of evaluations:  300  -2lnl: 1239.319519
 Number of evaluations:  400  -2lnl: 1239.303701
 Number of evaluations:  500  -2lnl: 1239.350357
 Number of evaluations:  600  -2lnl: 1239.417367
 Number of evaluations:  700  -2lnl: 1243.575689
 Number of evaluations:  800  -2lnl: 1239.311438
 Number of evaluations:  900  -2lnl: 1241.176613
 Number of evaluations:  1000  -2lnl: 1239.329909
 Number of evaluations:  1100  -2lnl: 1239.303675
cjs_m4
## 
## crm Model Summary
## 
## Npar :  12
## -2lnL:  1239.304
## AIC  :  1263.304
## 
## Beta
##                   Estimate
## Phi.(Intercept)  1.5987971
## Phi.time2       -2.6225405
## Phi.time3       -1.8076489
## Phi.time4       -1.1463136
## Phi.time5       -1.3796043
## Phi.time6       -0.6965196
## p.(Intercept)    1.2003709
## p.time3          1.5729892
## p.time4          1.2288969
## p.time5          0.7978441
## p.time6          1.2524201
## p.time7         -0.2257742
predict(object = cjs_m4)
## $Phi
##   time occ  estimate
## 1    6   6 0.7114173
## 2    5   5 0.5545798
## 3    4   4 0.6112295
## 4    3   3 0.4479760
## 5    2   2 0.2642989
## 6    1   1 0.8318502
## 
## $p
##   time occ  estimate
## 1    7   7 0.7260348
## 2    6   6 0.9207653
## 3    5   5 0.8806095
## 4    4   4 0.9190321
## 5    3   3 0.9412192
## 6    2   2 0.7685908

Multiple models + model selection (AIC)

When testing multiple hypotheses, you can speed up model formulation for multiple models. Caution! Models should always be biologically meaningful, don’t just fit too many models without thinking about them.

fit.starlings.cjs.models <- function() {
  # Apparent survival (Phi) formula
  Phi.time <- list(formula = ~ time)
  Phi.sex <- list(formula = ~ sex)
  Phi.sex.time <-
    list(formula = ~ sex * time) # interaction of sex and time
  Phi.habitat.time <- list(formula = ~ habitat * time)
  Phi.dot <- list(formula = ~ 1) # constant survival
  
  # Detection probability (p) formula
  p.sex <- list(formula = ~ sex)  # differs between males and females
  p.time <-
    list(formula = ~ time)  # one discrete estimate of p per capture event
  p.habitat <- list(formula = ~ habitat)
  p.dot <- list(formula = ~ 1) # constant detection
  
  # Construct all combinations and put into one model table
  cml <-
    create.model.list(c("Phi", "p")) # makes all possibile combinations of those parameter formulas
  results <- crm.wrapper(
    cml,
    data = m2_proc,
    ddl = m2_design,
    external = FALSE,
    accumulate = FALSE,
    hessian = TRUE
  )
  return(results)
}

# Run function
starlings_cjs_models <- fit.starlings.cjs.models()
## Phi.dot.p.dot 
## Phi.dot.p.habitat 
## 
 Number of evaluations:  100  -2lnl: 1291.647716Phi.dot.p.sex 
## 
 Number of evaluations:  100  -2lnl: 1289.126943Phi.dot.p.time 
## 
 Number of evaluations:  100  -2lnl: 1286.513634
 Number of evaluations:  200  -2lnl: 1284.979279
 Number of evaluations:  300  -2lnl: 1284.814101
 Number of evaluations:  400  -2lnl: 1284.787142
 Number of evaluations:  500  -2lnl:   1284.8424
 Number of evaluations:  600  -2lnl: 1285.724993
 Number of evaluations:  700  -2lnl:   1284.7819
 Number of evaluations:  100  -2lnl:  1284.88038
 Number of evaluations:  200  -2lnl: 1284.787694Phi.habitat.time.p.dot 
## 
 Number of evaluations:  100  -2lnl:  1227.93774
 Number of evaluations:  200  -2lnl: 1227.288078
 Number of evaluations:  300  -2lnl: 1227.014115
 Number of evaluations:  400  -2lnl: 1226.986396
 Number of evaluations:  500  -2lnl:  1226.97569
 Number of evaluations:  600  -2lnl: 1227.007945
 Number of evaluations:  700  -2lnl: 1226.992285
 Number of evaluations:  800  -2lnl: 1227.142454
 Number of evaluations:  900  -2lnl: 1226.984833
 Number of evaluations:  1000  -2lnl: 1227.016204
 Number of evaluations:  1100  -2lnl: 1226.985726
 Number of evaluations:  1200  -2lnl: 1227.030263
 Number of evaluations:  1300  -2lnl:  1226.99256
 Number of evaluations:  1400  -2lnl:  1226.97546
 Number of evaluations:  100  -2lnl: 1227.911352
 Number of evaluations:  200  -2lnl: 1227.036511
 Number of evaluations:  300  -2lnl: 1227.646509
 Number of evaluations:  400  -2lnl: 1227.025231
 Number of evaluations:  500  -2lnl: 1227.074533
 Number of evaluations:  600  -2lnl: 1226.990534
 Number of evaluations:  700  -2lnl: 1228.152346Phi.habitat.time.p.habitat 
## 
 Number of evaluations:  100  -2lnl:  1227.79551
 Number of evaluations:  200  -2lnl: 1227.120457
 Number of evaluations:  300  -2lnl: 1226.906674
 Number of evaluations:  400  -2lnl: 1226.799329
 Number of evaluations:  500  -2lnl: 1226.769009
 Number of evaluations:  600  -2lnl: 1226.766748
 Number of evaluations:  700  -2lnl: 1226.767176
 Number of evaluations:  800  -2lnl: 1226.839261
 Number of evaluations:  900  -2lnl: 1226.769448
 Number of evaluations:  1000  -2lnl: 1226.849578
 Number of evaluations:  1100  -2lnl: 1226.768915
 Number of evaluations:  1200  -2lnl: 1226.971117
 Number of evaluations:  1300  -2lnl: 1226.793747
 Number of evaluations:  1400  -2lnl:  1226.80454
 Number of evaluations:  1500  -2lnl: 1226.766644
 Number of evaluations:  100  -2lnl: 1227.604133
 Number of evaluations:  200  -2lnl: 1226.815512
 Number of evaluations:  300  -2lnl:  1227.74426
 Number of evaluations:  400  -2lnl:  1226.81178
 Number of evaluations:  500  -2lnl: 1226.980183
 Number of evaluations:  600  -2lnl: 1226.777627
 Number of evaluations:  700  -2lnl: 1228.149466
 Number of evaluations:  800  -2lnl: 1226.800166Phi.habitat.time.p.sex 
## 
 Number of evaluations:  100  -2lnl:  1226.29305
 Number of evaluations:  200  -2lnl: 1225.546506
 Number of evaluations:  300  -2lnl: 1225.331577
 Number of evaluations:  400  -2lnl: 1225.233333
 Number of evaluations:  500  -2lnl:  1225.22303
 Number of evaluations:  600  -2lnl: 1225.221478
 Number of evaluations:  700  -2lnl: 1225.222961
 Number of evaluations:  800  -2lnl: 1225.255474
 Number of evaluations:  900  -2lnl: 1225.242417
 Number of evaluations:  1000  -2lnl: 1225.280464
 Number of evaluations:  1100  -2lnl: 1225.239853
 Number of evaluations:  1200  -2lnl: 1225.252553
 Number of evaluations:  1300  -2lnl: 1225.232904
 Number of evaluations:  1400  -2lnl: 1225.239373
 Number of evaluations:  1500  -2lnl:  1225.22146
 Number of evaluations:  100  -2lnl: 1225.899994
 Number of evaluations:  200  -2lnl: 1225.267646
 Number of evaluations:  300  -2lnl: 1226.220945
 Number of evaluations:  400  -2lnl: 1225.270379
 Number of evaluations:  500  -2lnl: 1225.447483
 Number of evaluations:  600  -2lnl: 1225.232715
 Number of evaluations:  700  -2lnl: 1226.486483
 Number of evaluations:  800  -2lnl: 1225.258312Phi.habitat.time.p.time 
## 
 Number of evaluations:  100  -2lnl: 1236.170845
 Number of evaluations:  200  -2lnl: 1225.998227
 Number of evaluations:  300  -2lnl: 1224.510412
 Number of evaluations:  400  -2lnl: 1223.789346
 Number of evaluations:  500  -2lnl: 1223.614985
 Number of evaluations:  600  -2lnl: 1223.525763
 Number of evaluations:  700  -2lnl: 1223.408596
 Number of evaluations:  800  -2lnl: 1223.362659
 Number of evaluations:  900  -2lnl: 1223.314838
 Number of evaluations:  1000  -2lnl: 1223.310042
 Number of evaluations:  1100  -2lnl: 1223.333229
 Number of evaluations:  1200  -2lnl: 1224.496618
 Number of evaluations:  1300  -2lnl: 1223.314412
 Number of evaluations:  1400  -2lnl: 1223.694594
 Number of evaluations:  1500  -2lnl: 1223.327273
 Number of evaluations:  1600  -2lnl: 1223.897422
 Number of evaluations:  1700  -2lnl: 1223.326004
 Number of evaluations:  1800  -2lnl: 1223.548027
 Number of evaluations:  1900  -2lnl: 1223.314272
 Number of evaluations:  2000  -2lnl: 1223.381878
 Number of evaluations:  2100  -2lnl:  1223.32439
 Number of evaluations:  2200  -2lnl: 1223.714295
 Number of evaluations:  2300  -2lnl: 1223.314922
 Number of evaluations:  2400  -2lnl: 1223.331151
 Number of evaluations:  2500  -2lnl: 1223.310013
 Number of evaluations:  100  -2lnl: 1223.738771
 Number of evaluations:  200  -2lnl:  1223.39916
 Number of evaluations:  300  -2lnl: 1223.366651
 Number of evaluations:  400  -2lnl:  1223.43491
 Number of evaluations:  500  -2lnl: 1225.126663
 Number of evaluations:  600  -2lnl: 1223.602568
 Number of evaluations:  700  -2lnl: 1224.073011
 Number of evaluations:  800  -2lnl: 1223.329888
 Number of evaluations:  900  -2lnl: 1223.726268
 Number of evaluations:  1000  -2lnl: 1223.552908
 Number of evaluations:  1100  -2lnl: 1223.425592
 Number of evaluations:  1200  -2lnl:  1223.32932
 Number of evaluations:  1300  -2lnl:  1224.28456Phi.sex.p.dot 
## 
 Number of evaluations:  100  -2lnl: 1284.628321Phi.sex.p.habitat 
## 
 Number of evaluations:  100  -2lnl: 1284.362121
 Number of evaluations:  200  -2lnl: 1284.362121Phi.sex.p.sex 
## 
 Number of evaluations:  100  -2lnl: 1283.680512
 Number of evaluations:  200  -2lnl: 1283.680512Phi.sex.p.time 
## 
 Number of evaluations:  100  -2lnl: 1285.058065
 Number of evaluations:  200  -2lnl: 1278.590104
 Number of evaluations:  300  -2lnl: 1278.287167
 Number of evaluations:  400  -2lnl: 1278.550077
 Number of evaluations:  500  -2lnl: 1278.274914
 Number of evaluations:  600  -2lnl: 1278.278464
 Number of evaluations:  100  -2lnl: 1278.487128
 Number of evaluations:  200  -2lnl: 1278.281974Phi.sex.time.p.dot 
## 
 Number of evaluations:  100  -2lnl: 1226.568152
 Number of evaluations:  200  -2lnl: 1225.990052
 Number of evaluations:  300  -2lnl: 1225.526691
 Number of evaluations:  400  -2lnl: 1225.520238
 Number of evaluations:  500  -2lnl: 1225.441401
 Number of evaluations:  600  -2lnl: 1225.457658
 Number of evaluations:  700  -2lnl: 1227.229984
 Number of evaluations:  800  -2lnl: 1225.823925
 Number of evaluations:  900  -2lnl: 1226.666332
 Number of evaluations:  1000  -2lnl: 1225.557911
 Number of evaluations:  1100  -2lnl:  1226.41834
 Number of evaluations:  1200  -2lnl: 1225.546468
 Number of evaluations:  1300  -2lnl: 1225.441384
 Number of evaluations:  100  -2lnl: 1226.384775
 Number of evaluations:  200  -2lnl: 1225.590272
 Number of evaluations:  300  -2lnl: 1227.383442
 Number of evaluations:  400  -2lnl: 1225.783639
 Number of evaluations:  500  -2lnl:  1226.21125
 Number of evaluations:  600  -2lnl: 1225.580664
 Number of evaluations:  700  -2lnl: 1226.781155Phi.sex.time.p.habitat 
## 
 Number of evaluations:  100  -2lnl:  1226.27436
 Number of evaluations:  200  -2lnl: 1225.705522
 Number of evaluations:  300  -2lnl: 1225.343191
 Number of evaluations:  400  -2lnl: 1225.095631
 Number of evaluations:  500  -2lnl: 1225.060522
 Number of evaluations:  600  -2lnl: 1225.052203
 Number of evaluations:  700  -2lnl: 1225.070844
 Number of evaluations:  800  -2lnl: 1226.090349
 Number of evaluations:  900  -2lnl: 1225.070343
 Number of evaluations:  1000  -2lnl: 1225.434673
 Number of evaluations:  1100  -2lnl: 1225.075408
 Number of evaluations:  1200  -2lnl: 1225.277906
 Number of evaluations:  1300  -2lnl: 1225.063087
 Number of evaluations:  1400  -2lnl: 1225.427794
 Number of evaluations:  1500  -2lnl: 1225.059431
 Number of evaluations:  1600  -2lnl:  1225.05194
 Number of evaluations:  100  -2lnl: 1225.877588
 Number of evaluations:  200  -2lnl: 1225.472542
 Number of evaluations:  300  -2lnl: 1229.615687
 Number of evaluations:  400  -2lnl: 1225.117012
 Number of evaluations:  500  -2lnl: 1226.496245
 Number of evaluations:  600  -2lnl: 1225.122785
 Number of evaluations:  700  -2lnl: 1226.623784
 Number of evaluations:  800  -2lnl: 1225.094514Phi.sex.time.p.sex 
## 
 Number of evaluations:  100  -2lnl: 1226.094895
 Number of evaluations:  200  -2lnl:  1225.21013
 Number of evaluations:  300  -2lnl: 1224.774306
 Number of evaluations:  400  -2lnl: 1224.607338
 Number of evaluations:  500  -2lnl: 1224.548274
 Number of evaluations:  600  -2lnl: 1224.538384
 Number of evaluations:  700  -2lnl: 1224.787892
 Number of evaluations:  800  -2lnl: 1224.598502
 Number of evaluations:  900  -2lnl: 1225.403255
 Number of evaluations:  1000  -2lnl: 1224.579718
 Number of evaluations:  1100  -2lnl: 1224.871416
 Number of evaluations:  1200  -2lnl: 1224.548901
 Number of evaluations:  1300  -2lnl: 1225.183338
 Number of evaluations:  1400  -2lnl: 1224.563466
 Number of evaluations:  1500  -2lnl: 1224.538308
 Number of evaluations:  1600  -2lnl: 1224.538307
 Number of evaluations:  100  -2lnl: 1225.284145
 Number of evaluations:  200  -2lnl: 1225.049219
 Number of evaluations:  300  -2lnl: 1229.633306
 Number of evaluations:  400  -2lnl: 1224.615562
 Number of evaluations:  500  -2lnl: 1226.253282
 Number of evaluations:  600  -2lnl:  1224.62377
 Number of evaluations:  700  -2lnl: 1226.119317
 Number of evaluations:  800  -2lnl: 1224.588959Phi.sex.time.p.time 
## 
 Number of evaluations:  100  -2lnl: 1234.690315
 Number of evaluations:  200  -2lnl:  1225.62689
 Number of evaluations:  300  -2lnl: 1224.314752
 Number of evaluations:  400  -2lnl: 1223.075053
 Number of evaluations:  500  -2lnl: 1222.724628
 Number of evaluations:  600  -2lnl: 1222.508521
 Number of evaluations:  700  -2lnl: 1222.332531
 Number of evaluations:  800  -2lnl: 1222.210487
 Number of evaluations:  900  -2lnl: 1222.133208
 Number of evaluations:  1000  -2lnl: 1222.130014
 Number of evaluations:  1100  -2lnl: 1222.199949
 Number of evaluations:  1200  -2lnl: 1246.351238
 Number of evaluations:  1300  -2lnl: 1223.169168
 Number of evaluations:  1400  -2lnl: 1239.600774
 Number of evaluations:  1500  -2lnl: 1222.436111
 Number of evaluations:  1600  -2lnl: 1224.264099
 Number of evaluations:  1700  -2lnl: 1222.224401
 Number of evaluations:  1800  -2lnl:  1235.88525
 Number of evaluations:  1900  -2lnl: 1222.162762
 Number of evaluations:  2000  -2lnl: 1222.701721
 Number of evaluations:  2100  -2lnl: 1222.194462
 Number of evaluations:  2200  -2lnl: 1223.310846
 Number of evaluations:  2300  -2lnl:  1222.26185
 Number of evaluations:  2400  -2lnl: 1222.129981
 Number of evaluations:  2500  -2lnl: 1222.129981
 Number of evaluations:  100  -2lnl: 1222.614653
 Number of evaluations:  200  -2lnl: 1223.472304
 Number of evaluations:  300  -2lnl: 1225.226947
 Number of evaluations:  400  -2lnl: 1222.614453
 Number of evaluations:  500  -2lnl: 1224.253244
 Number of evaluations:  600  -2lnl: 1222.453988
 Number of evaluations:  700  -2lnl: 1226.374387
 Number of evaluations:  800  -2lnl: 1222.358573
 Number of evaluations:  900  -2lnl: 1225.989839
 Number of evaluations:  1000  -2lnl: 1222.382887
 Number of evaluations:  1100  -2lnl: 1222.225911
 Number of evaluations:  1200  -2lnl: 1222.161356
 Number of evaluations:  1300  -2lnl: 1223.262178Phi.time.p.dot 
## 
 Number of evaluations:  100  -2lnl: 1242.648575
 Number of evaluations:  200  -2lnl: 1242.487127
 Number of evaluations:  300  -2lnl: 1245.445377
 Number of evaluations:  400  -2lnl: 1242.418534
 Number of evaluations:  100  -2lnl: 1243.038192
 Number of evaluations:  200  -2lnl: 1242.446706Phi.time.p.habitat 
## 
 Number of evaluations:  100  -2lnl: 1242.360719
 Number of evaluations:  200  -2lnl: 1241.978159
 Number of evaluations:  300  -2lnl: 1243.089518
 Number of evaluations:  400  -2lnl: 1241.989515
 Number of evaluations:  500  -2lnl: 1242.041195
 Number of evaluations:  100  -2lnl:  1243.15084
 Number of evaluations:  200  -2lnl: 1242.099554Phi.time.p.sex 
## 
 Number of evaluations:  100  -2lnl: 1240.930825
 Number of evaluations:  200  -2lnl: 1240.477856
 Number of evaluations:  300  -2lnl: 1242.083998
 Number of evaluations:  400  -2lnl: 1240.508214
 Number of evaluations:  500  -2lnl: 1241.395697
 Number of evaluations:  100  -2lnl: 1241.640046
 Number of evaluations:  200  -2lnl: 1240.590849Phi.time.p.time 
## 
 Number of evaluations:  100  -2lnl: 1244.670707
 Number of evaluations:  200  -2lnl: 1239.725494
 Number of evaluations:  300  -2lnl: 1239.319519
 Number of evaluations:  400  -2lnl: 1239.303701
 Number of evaluations:  500  -2lnl: 1239.350357
 Number of evaluations:  600  -2lnl: 1239.417367
 Number of evaluations:  700  -2lnl: 1243.575689
 Number of evaluations:  800  -2lnl: 1239.311438
 Number of evaluations:  900  -2lnl: 1241.176613
 Number of evaluations:  1000  -2lnl: 1239.329909
 Number of evaluations:  1100  -2lnl: 1239.303675
 Number of evaluations:  100  -2lnl: 1247.504101
 Number of evaluations:  200  -2lnl: 1239.375924
 Number of evaluations:  300  -2lnl: 1240.074042
 Number of evaluations:  400  -2lnl: 1239.422373
 Number of evaluations:  500  -2lnl: 1239.433946
 Number of evaluations:  600  -2lnl: 1239.305696

Look at model results

summary(starlings_cjs_models)
##                            Length Class      Mode
## Phi.dot.p.dot              5      crm        list
## Phi.dot.p.habitat          5      crm        list
## Phi.dot.p.sex              5      crm        list
## Phi.dot.p.time             5      crm        list
## Phi.habitat.time.p.dot     5      crm        list
## Phi.habitat.time.p.habitat 5      crm        list
## Phi.habitat.time.p.sex     5      crm        list
## Phi.habitat.time.p.time    5      crm        list
## Phi.sex.p.dot              5      crm        list
## Phi.sex.p.habitat          5      crm        list
## Phi.sex.p.sex              5      crm        list
## Phi.sex.p.time             5      crm        list
## Phi.sex.time.p.dot         5      crm        list
## Phi.sex.time.p.habitat     5      crm        list
## Phi.sex.time.p.sex         5      crm        list
## Phi.sex.time.p.time        5      crm        list
## Phi.time.p.dot             5      crm        list
## Phi.time.p.habitat         5      crm        list
## Phi.time.p.sex             5      crm        list
## Phi.time.p.time            5      crm        list
## model.table                7      data.frame list

Time-dependent variables

Heat_wave <- matrix(rep(c(0, 1, 1, 0, 0, 0), each = nrow(starlings)), ncol =
                      6)
colnames(Heat_wave) = paste("Heat_wave", 1:6, sep = "")

starlings_heatwave <- cbind(starlings, Heat_wave)

#generate some random weight
starlings_heatwave$weight <- round(rnorm(nrow(starlings), 80, 3))




#prepare data
starlings.proc <- process.data(starlings_heatwave)

design.Phi <-
  list(static = c("weight", "sex"),
       time.varying = c("Heat_wave"))

design.p <- list(static = c("sex"))

design.parameters <- list(Phi = design.Phi, p = design.p)


m3_design <- make.design.data(data = starlings.proc,
                              parameters = design.parameters)

names(m3_design$Phi)
##  [1] "id"        "occ"       "time"      "cohort"    "age"       "Heat_wave"
##  [7] "weight"    "sex"       "Time"      "Cohort"    "Age"       "order"
names(m3_design$p)
##  [1] "id"     "occ"    "time"   "cohort" "age"    "sex"    "Time"   "Cohort"
##  [9] "Age"    "fix"    "order"
Phi.heat.sex <- list(formula =  ~ Heat_wave * sex)
p.sex <- list(formula =  ~ sex)

m3 <- crm(
  starlings.proc,
  m3_design,
  hessian = TRUE,
  model.parameters = list(Phi = Phi.heat.sex,
                          p = p.sex)
)
## 
 Number of evaluations:  100  -2lnl: 1244.121204
 Number of evaluations:  200  -2lnl: 1243.982328
 Number of evaluations:  300  -2lnl: 1244.263347
 Number of evaluations:  100  -2lnl: 1245.196681
predict(m3)
## $Phi
##   Heat_wave    sex occ  estimate         se       lcl       ucl
## 1         0 Female   6 0.5888999 0.03365102 0.5217320 0.6529102
## 2         0   Male   6 0.6029780 0.03050136 0.5419426 0.6609683
## 3         1 Female   3 0.2874165 0.03769815 0.2194742 0.3665152
## 4         1   Male   3 0.4682230 0.04526839 0.3813970 0.5570169
## 
## $p
##      sex occ  estimate         se       lcl       ucl
## 1 Female   7 0.8722827 0.03817459 0.7772436 0.9304042
## 2   Male   7 0.9191637 0.02691341 0.8482694 0.9585520

More complex models & predicting

Phi.weight.heat <-list(formula=~weight*Heat_wave)
p.sex <-list(formula=~sex)


m4 <- crm(
  starlings.proc,
  m3_design,
  hessian = TRUE,
  model.parameters = list(Phi = Phi.weight.heat,
                          p = p.sex)
)
## 
 Number of evaluations:  100  -2lnl: 1252.677436
 Number of evaluations:  200  -2lnl:  1255.89512
 Number of evaluations:  300  -2lnl: 1252.968609
 Number of evaluations:  400  -2lnl: 1252.676557
 Number of evaluations:  100  -2lnl: 1257.437678
#predict model output
predict(m4)
## $Phi
##    weight Heat_wave occ  estimate         se       lcl       ucl
## 1      79         0   6 0.6003154 0.02366211 0.5531702 0.6456726
## 2      78         0   6 0.6030211 0.02649719 0.5501114 0.6536262
## 3      76         0   6 0.6084135 0.03577190 0.5365250 0.6758874
## 4      87         0   6 0.5784621 0.05581375 0.4670070 0.6824580
## 5      84         0   6 0.5866981 0.03700753 0.5127970 0.6568898
## 6      77         0   6 0.6057205 0.03070027 0.5442351 0.6640317
## 7      82         0   6 0.5921622 0.02713700 0.5380946 0.6440871
## 8      80         0   6 0.5976036 0.02274455 0.5523348 0.6412679
## 9      81         0   6 0.5948858 0.02399419 0.5471243 0.6409161
## 10     85         0   6 0.5839579 0.04294503 0.4981462 0.6649651
## 11     75         0   6 0.6110998 0.04137545 0.5276484 0.6885113
## 12     83         0   6 0.5894329 0.03163288 0.5263314 0.6497225
## 13     86         0   6 0.5812125 0.04925269 0.4827940 0.6735658
## 14     73         0   6 0.6164522 0.05343803 0.5078826 0.7145325
## 15     74         0   4 0.6137795 0.04730548 0.5180168 0.7014801
## 16     90         0   4 0.5701823 0.07641024 0.4186053 0.7096513
## 17     80         1   3 0.3691765 0.02984352 0.3128306 0.4293305
## 18     86         1   3 0.4125508 0.06882940 0.2869833 0.5506308
## 19     78         1   3 0.3551391 0.03739509 0.2856581 0.4313161
## 20     77         1   3 0.3482111 0.04407812 0.2674567 0.4387446
## 21     81         1   3 0.3762808 0.03081215 0.3180633 0.4383054
## 22     79         1   3 0.3621284 0.03230885 0.3014750 0.4275162
## 23     84         1   3 0.3979045 0.04992160 0.3052082 0.4985524
## 24     83         1   3 0.3906476 0.04173233 0.3125535 0.4747780
## 25     82         1   3 0.3834388 0.03507171 0.3174010 0.4540764
## 26     76         1   3 0.3413467 0.05162750 0.2483609 0.4483786
## 27     87         1   3 0.4199341 0.07901220 0.2771182 0.5775476
## 28     85         1   3 0.4052065 0.05906622 0.2964706 0.5241130
## 29     75         1   2 0.3345482 0.05961159 0.2292578 0.4593738
## 30     88         0   6 0.5757068 0.06255564 0.4509656 0.6914950
## 31     72         0   5 0.6191180 0.05969651 0.4974060 0.7275026
## 32     74         1   3 0.3278177 0.06778289 0.2106540 0.4712427
## 33     72         1   2 0.3145687 0.08415756 0.1759748 0.4965423
## 
## $p
##      sex occ  estimate         se       lcl       ucl
## 1 Female   7 0.8568973 0.03983628 0.7600623 0.9188251
## 2   Male   7 0.9240681 0.02511068 0.8578462 0.9608486
new_starling <- expand.grid(sex=c("Female","Male"),weight=60:80,Heat_wave1=0,
       Heat_wave2=1,Heat_wave3=1,Heat_wave4=0,Heat_wave5=0,Heat_wave6=0) 


prediction <- predict(m4,newdata=new_starling,se=TRUE)
prediction$Phi$Heat_wave = factor(prediction$Phi$Heat_wave, labels = c("No_Heatwave", "Heatwave"))

ggplot(prediction$Phi, aes(weight, estimate, ymin = lcl, ymax = ucl)) +
  geom_errorbar(width = 0.2) +
  geom_point() +
  geom_line() +
  xlab("\nWeight") + ylab("Survival\n") + facet_grid(Heat_wave ~ .)

Exercise

Some exercises:

Exercise 6.1

The data set starlings contains one column (“infected”) that represents an infection with an incurable disease (i.e. a static variable comparable to sex). Formulate hypothesis associated with this infection status (write them down!) & try to build simple CMR models. Think about why infection could matter for capture probability or survival. There are some NAs in the data set, maybe try ?na-omit() and remove some of the observations without data from all analyses.Discuss the results in one sentence

Exercise 6.2

In our previous model (m3), we saw that sex & heatwave played an important role in individual survival. Formulate a model that additionally accounts for the habitat type (hint: * or :)

Exercise 6.3 (advanced)

This exercise will require some data-wrangling, the data set is not yet prepared: In the file locations, you will find measured temperatures for each year.

(Hint: You might need dplyr::left_join() and dplyr::rename() to get all data in the right format and bind the rows according to the habitat type.)

Test how temperature affects survival, and how sex affects detection probability. Then predict the survival probabilities for extreme weather scenarios (temperatures of 30,33,35,38,41,44 degree Celsius)

Solution 6.1

There are many solutions, one is:

Hypothesis: infection affects survival & prediction

Prediction: Lower survival AND detection probability in infected individuals. Here, we hypothesize that infected animals move less, therefore are less likely to be captured. Additionally, infected animals suffer from disease, thus survival is decreased

package.list=c("here",
               "marked",
               "tidyverse")

for (package in package.list) {
  if (!require(package, character.only=T, quietly=T)) {
    install.packages(package)
    library(package, character.only=T)
  }
}


starlings <-
  read_csv(here("data", "data_berlin", "animal_data", "starlings.csv")) %>%
  mutate(
    sex = as.factor(sex),
    habitat = as.factor(habitat),
    infected = as.factor(infected)
  ) %>%
  rename(ID = "...1")%>%
  as.data.frame()

#load environmental data on capturing
locations <-
  read_csv(here(
    "data",
    "data_berlin",
    "animal_data",
    "starlings_capture_location_data.csv"
  ))


temperatures <- locations %>%
  dplyr::select(-c("Lon", "Lat")) %>%
  rename(habitat = Habitat)


m_infection<- process.data(data = na.omit(starlings),groups = c("infected"))
m_infetion_design <- make.design.data(data = m_infection)

phi_infection <- list(formula = ~ infected)
p_infection <- list(formula = ~ infected)

cjs_infection <- crm(
  m_infection,
  m_infetion_design,
  model.parameters = list(Phi = phi_infection,
                          p = p_infection),
  accumulate = FALSE
)
## 
 Number of evaluations:  100  -2lnl: 861.9151509
cjs_infection
## 
## crm Model Summary
## 
## Npar :  4
## -2lnL:  861.9147
## AIC  :  869.9147
## 
## Beta
##                    Estimate
## Phi.(Intercept) -0.03261002
## Phi.infectedYes  0.20133675
## p.(Intercept)    2.68053887
## p.infectedYes   -0.76975772
predict(object = cjs_infection)
## $Phi
##   infected occ  estimate
## 1      Yes   6 0.5420819
## 2       No   6 0.4918482
## 
## $p
##   infected occ  estimate
## 1      Yes   7 0.8711069
## 2       No   7 0.9358685

Solution 6.2

#prepare data
Heat_wave <- matrix(rep(c(0, 1, 1, 0, 0, 0), each = nrow(starlings)), ncol =
                      6)
colnames(Heat_wave) = paste("Heat_wave", 1:6, sep = "")

starlings_heatwave <- cbind(starlings, Heat_wave)

#generate some random weight
starlings_heatwave$weight <- round(rnorm(nrow(starlings), 80, 3))


starlings.proc <- process.data(starlings_heatwave)

design.Phi <-
  list(static = c("weight", "sex","habitat"),
       time.varying = c("Heat_wave"))

design.p <- list(static = c("sex"))

design.parameters <- list(Phi = design.Phi, p = design.p)


m6_design <- make.design.data(data = starlings.proc,
                              parameters = design.parameters)

names(m6_design$Phi)
##  [1] "id"        "occ"       "time"      "cohort"    "age"       "Heat_wave"
##  [7] "weight"    "sex"       "habitat"   "Time"      "Cohort"    "Age"      
## [13] "order"
names(m6_design$p)
##  [1] "id"     "occ"    "time"   "cohort" "age"    "sex"    "Time"   "Cohort"
##  [9] "Age"    "fix"    "order"
Phi.heat.sex.habitat <- list(formula =  ~ Heat_wave * sex * habitat)
p.sex <- list(formula =  ~ sex)

m6 <- crm(
  starlings.proc,
  m6_design,
  hessian = TRUE,
  model.parameters = list(Phi = Phi.heat.sex.habitat,
                          p = p.sex)
)
## 
 Number of evaluations:  100  -2lnl:  1224.55332
 Number of evaluations:  200  -2lnl: 1224.114866
 Number of evaluations:  300  -2lnl:  1224.07174
 Number of evaluations:  400  -2lnl: 1224.076663
 Number of evaluations:  500  -2lnl: 1224.093993
 Number of evaluations:  600  -2lnl: 1224.305982
 Number of evaluations:  700  -2lnl: 1224.072706
 Number of evaluations:  800  -2lnl: 1224.071643
 Number of evaluations:  100  -2lnl: 1224.265117
 Number of evaluations:  200  -2lnl: 1224.137044
 Number of evaluations:  300  -2lnl: 1224.878393
 Number of evaluations:  400  -2lnl: 1224.083587
predict(m6)
## $Phi
##   Heat_wave    sex habitat occ  estimate         se        lcl       ucl
## 1         0 Female   rural   6 0.5359784 0.04604001 0.44554937 0.6241018
## 2         0   Male   rural   6 0.5956507 0.04259539 0.51018411 0.6756843
## 3         1 Female   rural   3 0.4429930 0.05930527 0.33178993 0.5602196
## 4         1   Male   rural   3 0.4807007 0.06379440 0.35936302 0.6043584
## 5         0 Female   urban   6 0.6391793 0.04551414 0.54611911 0.7228425
## 6         0   Male   urban   6 0.6103285 0.04248132 0.52462849 0.6897163
## 7         1 Female   urban   3 0.1365968 0.04010417 0.07513682 0.2355269
## 8         1   Male   urban   3 0.4555735 0.06396440 0.33544774 0.5811012
## 
## $p
##      sex occ  estimate         se       lcl       ucl
## 1 Female   7 0.8739003 0.03775411 0.7797819 0.9313357
## 2   Male   7 0.9191464 0.02691835 0.8482405 0.9585424

Solution 6.3

starlings_temperatur <- starlings %>%
  left_join(temperatures, by = "habitat")%>%
  rename(temperature1=temperature_2016,
         temperature2=temperature_2017,
         temperature3=temperature_2018,
         temperature4=temperature_2019,
         temperature5=temperature_2020,
         temperature6=temperature_2021,
         temperature7=temperature_2022)

#prepare data
starlings.proc.temp <- process.data(starlings_temperatur)

design.Phi <-
  list(static = c("sex"),
       time.varying = c("temperature"))

design.p <- list(static = c("sex"))

design.parameters <- list(Phi = design.Phi, p = design.p)


m7_design <- make.design.data(data = starlings.proc.temp,
                              parameters = design.parameters)

names(m7_design$Phi)
##  [1] "id"          "occ"         "time"        "cohort"      "age"        
##  [6] "temperature" "sex"         "Time"        "Cohort"      "Age"        
## [11] "order"
names(m7_design$p)
##  [1] "id"     "occ"    "time"   "cohort" "age"    "sex"    "Time"   "Cohort"
##  [9] "Age"    "fix"    "order"
Phi.heat.sex <- list(formula =  ~ temperature)
p.sex <- list(formula =  ~ sex)

m7 <- crm(
  starlings.proc.temp,
  m7_design,
  hessian = TRUE,
  model.parameters = list(Phi = Phi.heat.sex,
                          p = p.sex)
)
## 
 Number of evaluations:  100  -2lnl: 1249.569476
 Number of evaluations:  200  -2lnl: 1242.783021
predict(m7)
## $Phi
##   temperature occ  estimate         se       lcl       ucl
## 1        29.0   6 0.5900881 0.02104224 0.5483047 0.6306103
## 2        27.0   4 0.6457556 0.02513997 0.5950862 0.6933517
## 3        31.0   3 0.5320125 0.01892171 0.4948173 0.5688552
## 4        35.0   2 0.4148422 0.02469059 0.3674108 0.4639063
## 5        29.5   1 0.5757385 0.02023851 0.5356662 0.6148391
## 6        28.0   5 0.6183101 0.02298216 0.5723670 0.6622301
## 7        30.0   4 0.5612600 0.01959489 0.5225627 0.5992263
## 8        36.0   3 0.3865030 0.02722904 0.3346766 0.4410349
## 9        43.0   2 0.2161236 0.03866961 0.1498497 0.3013197
## 
## $p
##      sex occ  estimate         se       lcl       ucl
## 1 Female   7 0.8699031 0.03750525 0.7773961 0.9275499
## 2   Male   7 0.9266009 0.02440500 0.8620348 0.9622732
new_starling2 <- expand.grid(sex=c("Female","Male"),temperature1=30,
       temperature2=33,temperature3=35,temperature4=38,temperature5=41,temperature6=44,weight=60:80) 

prediction <- predict(m7,newdata=new_starling2,se=TRUE)
prediction$Phi$temperature = factor(prediction$Phi$temperature)


ggplot(prediction$Phi, aes(temperature, estimate, ymin = lcl, ymax = ucl)) +
  geom_errorbar(width = 0.2) +
  geom_point() +
  geom_line() +
  xlab("\nTemperature") + ylab("Survival probability\n")

Session Info 
Sys.time()
## [1] "2023-03-16 19:09:26 CET"
#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=C                              
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] tmap_3.3-3      sf_1.0-9        forcats_0.5.2   stringr_1.5.0  
##  [5] dplyr_1.0.10    purrr_1.0.1     readr_2.1.3     tidyr_1.3.0    
##  [9] tibble_3.1.8    ggplot2_3.4.0   tidyverse_1.3.2 skimr_2.1.5    
## [13] marked_1.2.6    lme4_1.1-31     Matrix_1.5-3    here_1.0.1     
## 
## loaded via a namespace (and not attached):
##   [1] leafem_0.2.0            googledrive_2.0.0       minqa_1.2.5            
##   [4] colorspace_2.1-0        ellipsis_0.3.2          class_7.3-21           
##   [7] leaflet_2.1.1           rprojroot_2.0.3         base64enc_0.1-3        
##  [10] fs_1.6.0                dichromat_2.0-0.1       rstudioapi_0.14        
##  [13] proxy_0.4-27            farver_2.1.1            bit64_4.0.5            
##  [16] fansi_1.0.4             lubridate_1.9.1         xml2_1.3.3             
##  [19] codetools_0.2-18        splines_4.2.2           cachem_1.0.6           
##  [22] knitr_1.41              jsonlite_1.8.4          nloptr_2.0.3           
##  [25] tmaptools_3.1-1         broom_1.0.2             dbplyr_2.3.0           
##  [28] png_0.1-8               compiler_4.2.2          httr_1.4.4             
##  [31] backports_1.4.1         assertthat_0.2.1        fastmap_1.1.0          
##  [34] gargle_1.2.1            cli_3.6.0               s2_1.1.2               
##  [37] leaflet.providers_1.9.0 htmltools_0.5.4         tools_4.2.2            
##  [40] coda_0.19-4             gtable_0.3.1            glue_1.6.2             
##  [43] wk_0.7.1                Rcpp_1.0.10             raster_3.6-14          
##  [46] cellranger_1.1.0        jquerylib_0.1.4         vctrs_0.5.2            
##  [49] nlme_3.1-161            leafsync_0.1.0          crosstalk_1.2.0        
##  [52] lwgeom_0.2-11           xfun_0.36               rvest_1.0.3            
##  [55] timechange_0.2.0        lifecycle_1.0.3         XML_3.99-0.13          
##  [58] googlesheets4_1.0.1     terra_1.7-3             MASS_7.3-58.2          
##  [61] scales_1.2.1            vroom_1.6.1             ragg_1.2.5             
##  [64] hms_1.1.2               expm_0.999-7            TMB_1.9.2              
##  [67] RColorBrewer_1.1-3      yaml_2.3.7              sass_0.4.5             
##  [70] stringi_1.7.12          highr_0.10              e1071_1.7-12           
##  [73] optimx_2022-4.30        boot_1.3-28.1           truncnorm_1.0-8        
##  [76] R2admb_0.7.16.3         repr_1.1.5              rlang_1.0.6            
##  [79] pkgconfig_2.0.3         systemfonts_1.0.4       evaluate_0.20          
##  [82] lattice_0.20-45         labeling_0.4.2          htmlwidgets_1.6.1      
##  [85] bit_4.0.5               tidyselect_1.2.0        magrittr_2.0.3         
##  [88] bookdown_0.32           R6_2.5.1                generics_0.1.3         
##  [91] DBI_1.1.3               pillar_1.8.1            haven_2.5.1            
##  [94] withr_2.5.0             units_0.8-1             stars_0.6-0            
##  [97] sp_1.6-0                abind_1.4-5             modelr_0.1.10          
## [100] crayon_1.5.2            KernSmooth_2.23-20      utf8_1.2.2             
## [103] tzdb_0.3.0              rmarkdown_2.20          grid_4.2.2             
## [106] readxl_1.4.1            data.table_1.14.6       rmdformats_1.0.4       
## [109] reprex_2.0.2            digest_0.6.31           classInt_0.4-8         
## [112] numDeriv_2016.8-1.1     textshaping_0.3.6       munsell_0.5.0          
## [115] viridisLite_0.4.1       bslib_0.4.2